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Abstract 



In this work, we consider the problem of reconstructing a small anomaly 
in a viscoelastic medium from wave-field measurements. We choose Szabo's 
model to describe the viscoelastic properties of the medium. Expressing the 
ideal elastic field without any viscous effect in terms of the measured field in a 
^ . viscous medium, we generalize the imaging procedures, such as time reversal, 

Kirchhoff Imaging and Back propagation, for an ideal medium to detect an 
anomaly in a visco-elastic medium from wave-field measurements. 



1 Introduction 

>. 

^^ ^ We consider the problem of reconstructing a small anomaly in a viscoelastic medium 

0^ ■ from wave-field measurements. The Voigt model is a common model to describe the 

viscoelastic properties of tissues. Catheline et al. [TD] have shown that this model 
is well adapted to describe the viscoelastic response of tissues to low-frequency 

^^ [ excitations. We choose a more general model derived by Szabo et al. [TB] that 

describes observed power-law behavior of many viscoelastic materials. It is based 
on a time-domain statement of causality [15] . It reduces to the Voigt model for the 
specific case of quadratic frequency loss. Expressing the ideal elastic field without 
any viscous effect in terms of the measured field in a viscous medium, we generalize 

^ ■ the methods described in [H [31 HJ [SJ [5] ; namely the time reversal, back-propagation 

JH I and Krichhoff Imaging, to recover the viscoelastic and geometric properties of an 

anomaly from wave-field measurements. 

The article is organized as follows. In section [5] we introduce a general visco- 
elastic wave equation, section [3] is devoted to the derivation of the Green function 
in a viscoelastic medium. In section [5] we present anomaly imaging procedures and 
reconstruction methods in visco-elastic media. Numerical illustrations are provided 
in section [5l 

2 General Visco-Elastic Wave Equation 

When a wave travels through a biological medium, its amplitude decreases with 
time due to attenuation. The attenuation coefficient for biological tissue may be 
approximated by a power-law over a wide range of frequencies. Measured attenua- 
tion coefficients of soft tissue typically have linear or greater than linear dependence 
on frequency [11] [TSl [H] . 



*Centre de Mathematiqucs Appliquccs, CNRS UMR 7641, Ecolc Polytechniquc, 91128 
Palaiseau, France (bretin@cmap.polytechnique.fr, lili.guadarrama-bustos@cmap.polytechnique.fr, 
wahab@cmap.polytechnique.fr). 



In an ideal medium; without attenuation, Hooke's law gives the following rela- 
tionship between stress and strain tensors: 

T = C:S (1) 

where T, C and S are respectively stress, stiffness and strain tensors of orders 2, 4 
and 2 and : represents tensorial product. 

Consider a dissipative medium. Suppose that the medium is homogeneous and 
isotropic. We write 

C = [Cijki] = [^6ij5ki + fj-{5ik5ji + SiiSjk)] , (2) 

V = [Vijki] = blsSijSki + ripiSikSji + SiiS-jk)] , (3) 

where Sat is the Kronccker delta function, ^, A are the Lame parameters, and 77^, rjp 
are the shear and bulk viscosities, respectively. Here we have adopted the general- 
ized summation convention over the repeated index. 
Throughout this work we suppose that 

Vp,Vs«l- (4) 

For a medium obeying a power-law attenuation model and under the smallness 
condition ^, a generalized Hooke's law reads [TB] 

Tix,t)=C ■.S{x,t)+Tj:M{S)ix,t) (5) 

where the convolution operator A4 is given by 

_(_l)j//2^__5 y jg g^jj^ even integer, 

M{S)= < f (y - l)!(-l)fe+i)/2£a ^s y is an odd integer, (6) 

-|r{2/) sin(y7r/2)-^ * 5 y is a non integer. 

Here H{t) is the Heaviside function and T denotes the gamma function. 

Note that for the common case, y = 2, the generalized Hooke's law (O reduces 
to the Voigt model, 

r = C:5 + ,:-. (7) 

Taking the divergence of ([5]) we get 

V • r = (A + /2) V(V • u) + ^Au, 
where 

X = X + r]pM{-) and p. = fi + t]sM{-). 
Next, considering the equation of motion for the system, i.e., 

with p being the constant density and F the applied force. Using the expression for 
V • T, we obtain the generalized visco-elastic wave equation 

P-g^-F={\ + fl)ViV-u) + flAu. (9) 



3 Green's Function 

In this section we find the Green function of the viscoelastic wave equation ©. For 
doing so, we first need a Helmholtz decomposition. 

3.1 Helmholtz Decomposition 

The following lemma holds. 

Lemma 3.1 // the displacement field u{x,t) satisfies (^, g^' ' ~ VA + V x i? 
and u{x, 0) = VC + V x D and if the body force F = V<p/ +\^ xtpj then there exist 
potentials ip^ and ip^^ such that 

• u = Vifu + V X Vu,- V • -0M = 0; 



^ = f + c^A^„ + ^sM{A^.) « ^ - ^^^^ + clA^u + ^A^(a,V«), 



p^t 



with 



\ + 2^ 



2 P 



Vp + 27?s 



P P 

Proof. FoT ipu and -0„ defined as 



Vuix,t) 



i'u{x,t) 



^0 

t I'T 



-^ + (c2 + j/pAl)(V-u) 



1/;/ 



(c^ + i/,,7W)(V X u) 



and Us ~ — . 
P 




dsdr + tA + C 


(10) 


dsdr + tB + D 


(11) 



we have the required expression for u. Moreover, it is evident from (|lip that 
V • ^„ == 

Now, on differentiating ipu and ip^ twice with respect to time, we get 

"^^ - '^^ + 4A^u + UpM{A^u) 



df^ P 

Ot'' p 

Finally, applying Ai on last two equations, neglecting the higher order terms in i/g 
and Up and injecting back the expressions for M{Aipu) and A^(A7/'„), we get the 
required differential equations for tpu and tpu- n 



Let 



Kmii^) ^ uj\ [I 5-A^(a;) , m = s,p. 



(12) 



where the multiplication operator A4{uj) is the Fourier transform of the convolution 
operator M. 

If ifu and ^jju are causal then it implies the causality of the inverse Fourier 
transform of Km{^),m = s,p. Applying the Kramers-Kronig relations, it follows 
that 



^mKrn{(^) — T~L 



5ReX™(c^) 



and '^eK„i[ijj) = H 



'^rnKmi,'^) 



, m = p, s, 



(13) 



where T-L is the Hilbert transform. Note that Ti,"^ = —I. The convohition operator 
Ai given by (JH) is based on the constraint that causahty imposes on ([5]) . Under the 
smallness assumption (j3]), the expressions in ^ can be found from the Kraniers- 
Kronig relations ([T^ . One drawback of P^ is that the attenuation, '^mKm{u>), 
must be known at aU frequencies to determine the dispersion, 3fie_ftr,„(a;). However, 
bounds on the dispersion can be obtained from measurements of the attenuation 
over a finite frequency range |13j . 

3.2 Solution of ([9]) with a Concentrated Force. 

Let Uij denote the i-th component of the solution u^ of the elastic wave equation 
related to a force F concentrated in the Xj-direction. Let j = 1 for simplicity and 
suppose that 

F = -Tit)d{x - Oei - -T{t)d{x - 0(1, 0, 0), (14) 

where ^ is the source point and (61,62,63) is an orthonormal basis of R'^. The 
corresponding Helmholtz decomposition of the force F can be written [Tl] as 

' F = Vipf + V X -0/, 
V'/ = ^alr(7)' (15) 

^f ~ ~~A^ y^ d^ w) '~a57 w/J ' 

where r = |a; — ^|. 

Consider the Helmholtz decomposition for u^i as 

Mil = VV31 + V X -01 (16) 

where ipi and tpi are the solutions of the equations 



A^.-i^ + SA^c?*.)-^^^^^-*- m 



Cp ot- c^ pel c^p 



^^1 - -^ + -,M{d-t^i) = ^-^^^^^ ~ ^. (18) 



Taking the Fourier transform of ([T5| . ([T7)) and p^ with respect to t we get 

ui = V(^i + V X 01 (19) 

A^. + ^^. ^ "^^^^^^ - % (20) 

A^, + ^i.,= '^f^^-% (21) 

c; pet pcj 

with i^„i(w), m = p,s, given by p^ . 

It is well known that the Green's functions of the Helmholtz equations ([^D|) and 
(PIj) are 

5™(r,w) = , m = s,p. 

Awr 

Thus, following [Tl] we write (pi as 

\ S / P(47rcp)^ Jy dxi \x-Q 



and divide the volume V into spherical shells of radius h centered at observation 
point X. On each shell g^{x — Xi w) rests constant. So we have 

with h = \x — x\i R — lx~?l ^-nd da the appropriate surface element. 

Asm 

f d f l\ fO h>r 



Therefore, we have following expression for ipi: 
In the same way, the vector ipi is given by 
Introduce the following notation: 



(23) 



I„,{x, oj) = A„, r "" Ce^^^'"^"^'^ dC (24) 

Jo 

EmXx,io) = Ame^^''"^^'^^^, (25) 

/ i^™7W(tj)\ 
An(w) = 1 2 ' 'm = p,s. (26) 

We obtain, after a lengthy but simple calculation, that un is given by 

and therefore, it follows that the solution Uij for an arbitrary j is 

"y = W (37»7j - %) ;ir [/s(r,a;) - /p(r,a;)] + ^^jaj^^pir,^) 

+ 1^ (% - 7*7j) 7^s(?','^), 
where 7^ = {xi - Ci)/r. 

3.3 Green's function 

If we substitute T{t) = S{t), where delta is the Dirac mass, then the function u.y = 
Gij is the i-th component of the Green function related to the force concentrated 
in the Xj -direction. In this case, we have T{ll!) = 1. Thus, we have the following 
expression for Gij : 

^■u == 4^ (37ai - Sij) -p, [Isir,uj) - Ipir,uj)] + j^JHj^ Ep{r,uj) 



and 



which imphcs that 

G,,(r,w;0 -^.'^■(^w) +5l,(r,w) +C('^^)' (27) 

where 

9T,{rM = ^ (37.7, - hi) ^ [Is{r,co) /,(r,c.)] , (28) 

gP^{r,u;) = ^^ajg'ir^u;), (29) 

9h(r,u^) = ^ {S., - 7.7.)r(^,^)- (30) 

Let G{r,t;£^) = {Gij{r,t;^)) denote the transient Green function of © associ- 
ated with the source point £,. Let G'^{r,t;£^) and Wm{r,t) be the inverse Fourier 
transforms of A„i(a;)5'"(r, w) and Im{r^Lo)^m = p,s, respectively. Then, from (j27t - 
[50)) . we have 

G,j(r, t; = ^1^1, G^r, i; + ^ (-5,, - 7.7,) G^(^ i; 

(31) 
+ 4^^ (37.7, - 5.j) ^ [Ws{r, t) - Wp{r, t)] . 

Note that by a change of variables, 

4^ r ^2^m, 



VF„(r,i) = — / C'G"(C,i;e)dC- 

4 Imaging procedure 

Consider the limiting case A — > +oo. The Green function for a quasi-incompressible 
visco-elastic medium is given by 

G.,(r,i;0 = ^(%-7.7,)G«(r-,t;0 

+^ (37.7, - hj) ^ Jo eC^iC t; OdC 

To generalize the detection algorithms presented in [21 [3l [3 IH |8] to the visco-elastic 
case we shall express the ideal Green function without any viscous effect in terms 
of the Green function in a viscous medium. From 

G-'(r,i;0 = -^ / e-^^'^'A,{u;)g'{r,uj)duj, 
V Ztt J a 



it follows that 

G%r,t;0 = ^= / Mu^) dio. 

V 27r Jm 47rr 



4.1 Approximation of the Green Function 

Introduce the operator 





1 p /•+00 



for a causal function 6. Wc have 



G^(M;e) = i.(fc^^), 



Anr 



and therefore, 



L*G''{r,t-Cj ^L*L{ 



5{t ~r/cs) . 



where L* is the L^{0, +oo)-adjoint of L. 

Consider for simphcity the Voigt modeL Then, Ai{uj) = — \/— Iw and hence, 



ll^s 



Ks{lo) = wW 1 H 5 — w ss w + 



-li/, 



s 2 



under the smahness assumption Q. The operator L can then be approximated by 



ZTT 



+ CK3 



A,(w)(/)(T)e" 



-i"(^-*)drda;. 



Since 



and 






/J/^T 






it follows that 



Analogously, 



L0(t) 



+ CXD 



dt 



t C CsC-r-t) 

' V) — P 2..3T (^7-_ 



r \Pjuj7r 






(32) 



(33) 



Since the phase in (j33]) is quadratic and Us is small then by the stationary phase 
theorem lA.li we can prove the following theorem: 



Theorem 4.1 



L*0«0+^au(i0), U^<l,^^tdu^ 



and therefore, 



L*L<pK(^+^dt{tdt<p), 



iL*Ly'(f>^(f>-^dtitdt<f>). 



(34) 
(35) 



Proof. (See appendix 1X1) 



4.2 Reconstruction Methods 

From the previous section, it follows that the ideal Green function, S{t — r / Cs) / {4:TTr) , 
can be approximately reconstructed from the viscous Green function, C^lr, t; ^), by 
either solving the ODE 

cb+^dt{tdtcl>)=L*G'{r,t;0, 

with = 0, t ^ or just making the approximation 

S{t - r/cs)/{4TTr) ^ L*G-'{r,t;0 ' '^dt{tdtL*G'{r,t;0). 

c 

Once the ideal Green function (5(r — r/cs)/(47rr) is reconstructed, one can find 
its source ^ using a time-reversal, a Kirchhoff or a back-propagation algorithm. See 

[iisiiiig. 

Using the asymptotic formalism developed in [SJ O [7], one can also find the 
shear modulus of the anomaly using the ideal near-field measurements which can 
be reconstructed from the near-field measurements in the viscous medium. The 
asymptotic formalism reduces the anomaly imaging problem to the detection of the 
location and the reconstruction of a certain polarizability tensor in the far-field and 
separates the scales in the near-field. 

5 Numerical Illustrations 

In this section, we illustrate the profile of the Green function. We choose parameters 
of simulation as in the work of Bercoff et al.^: we take p = 1000, Cg = 1, Cp = 40, 
r = 0.015 and Vp = 0. 



In figure 1, we plot temporal representation of the green function: 
pCp Anpr 



t ^ — (G^r^tiO + G'{r,t;0) + -. ^ [Ws{r,t) - Wp{r,t)] 



for three different values of y and Vs- We can see that the attenuation behavior 
varies with respect to different choices of power law exponent y. One can clearly 
distinguish the three different terms of the Green function; i. e. Gf, , Gf, and Gf ' . 
Figure 2 corresponds to spatial representation of the green function: 

(x,2/) ^ -^ {{xlrfGP{r,t-i) + (1 - (a;/r)2)G^(r,i; 0)+^(3(a;/r)2-l) [W.{r,t) - Wp{r,t)] 

for different values of y at t = 0.015. As expected, we get a diffusion of the wavefront 
with the increasing values of y and depending the choice of Vs . 

In figure 3, we illustrate the results of the approximation of the operator Lc/i 
with the smooth function 0(t) — exp(— 50 * (t — 1).^)". As shown by the stationary 
phase thcorcm lA.il . the numerically calculated L°°-error 



l|i^-(0+^¥')llL~(K+) 



is of order two. 
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Figure 1: Temporal response to a spatio-temporal delta function using a purely 
elastic Green's function (red line) and a viscous Green's function (blue line): Left, 
y = 1.5, I/, = 4 ; Center, y^2,v,^ 0.2 ; Right, y = 2.5, v,, = 0.002. 
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Figure 2: 2D spatial response to a spatio-temporal delta function at ^ = 0.015 with 
a purely elastie Green's function, a viscous Green's function with y = 2, z/g = 0.2 
and y = 2.5, Vs = 0.002. 
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Figure 3: Approximation of L via stationary phase theorem : Left, comparison 
between L4> and + ^txf/' where -jf = 0.0001 and is a smooth function. Right: 
error % -^ ||i0 — </> + %||oo in logarithmic scale. 

6 Conclusion 

In this paper, we have computed the Green function in a visco-elastic medium obey- 
ing a frequency power-law. For the Voigt model, which corresponds to a quadratic 
frequency loss, we have used the stationary phase theorem lA.ll to reconstruct the 
ideal Green function from the viscous one by solving an ODE. Once the ideal Green 
function is reconstructed, one can find its source ^ using the algorithms in [H |31 HI [S] 
such as time reversal, back-propagation, and Kirchhoff Imaging. For more general 
power-law media, one can recover the ideal Green function from the viscous one by 
inverting a fractional derivative operator. This would be the subject of a forthcom- 
ing paper. 



A Proof of Theorem (14.11) 



The proof of theorem (|4.ip is based on the following theorem (see [T^ Theorem 

7.7.1]). 

Theorem A.l (Stationary Phase)Lei K C [0, oo) he a compact set, X an open 
neighborhood of K and k a positive integer. If t/j G C'^''{K), f G C^^^^{X) and 
Im(f) >QinX, Im{f{to)) = 0, /'(io) = 0, /"(to) ^ 0, f ^ tn K \ {to} then 
for e> 



K 



^{t)e'f'^*^/'dx - e'f<^*°^^' (A/"(io)/27ri)~^^^ J2 ^^^J^ 

j<k 



<Ce^ Y^ sup|V'^"^(a;)|. 



a<2k 



Here C is hounded when f stays in a bounded set in C^'^'*"^(A") and \t — to\/\f'{t)\ 
has a uniform bound. With, 

9toit) = f{t)-.f{to)-lf"{to){t-to)\ 
which vanishes up to third order at to ; we have 

L^i'= E E *-^'^(-i)v"(to)-''«V')('^nio). □ 



10 



Note that Li can be expressed as the sum Liil' = L\il} + L\i}} + if V': where L\ is 
respectively associate to the pair (vj, fij) = (1, 0), (2, 1), (3, 2) and is identified to 

Now we turn to the proof of formula ([M|. Let us first consider the case of 
operator L* . We have 



L*m 



+ CXD 



-Mr) ""' 



t \j2jujsi 



1 / /•+°° 



e--2i77i-dT^—^{ I ^(r)e*^(^)/' 



with, /(r) = i7r(T — i)^, e = ^^a°* and V'(''') = t(J){t). Remark that the phase / 
satisfies at r = t , /(t) = 0, f'{t)'= 0, /"(t) = 2i7r 7^ 0. Moreover, we have 

'e»/(*)A(e-i/"(t)/2j^)-i/2 = y^ 
5t(r) = /(r)-/(t)-i/"W(r-t)2 = 

Thus, Theorem I A. 1 1 implies that 



The case of the operator L is very similar. Note that 



Uit) 



+00 






+00 






with/(T)=*7r(l^,e-^ 



2^^2 and V'(''') = 4>{t)t 2 . It follows that 
/'(r) = *7r ( 1 - ^ V /"(^)^2»7r^, /"(i) = 2»7ri 



and the function (7t(r) equals to 



, (r - t)^ . (r - t)^ . (i - rf 

gt[Tj = ITT in = ZTT- 



Tt 



We deduce that 






and then, 



i?^ = ^/"(r^ (fff w^w + 45^ wv^'w) = ^ (3 (^)' - 3|i) = ^ (3^ 

where 4'{t) = (J){t)/t. Then, we have 



9 0(t) 

2 (3/2 



11 



LV = Lli; + Lli; + LI^ 



^U"«.,.-.,^.(^^f)fiu_±^(«o)"^^m 



and again Thcoreni lA.il shows that 



L0(t) - { m + ^i4>"{t) 



<CTe3/2^sup|V''"Hi)l- 
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